clear all
version 16
cd "____"
set more off


// Some city/years with very good air quality may not converge with Stata automatic starting value. Changing starting value to (0.31,-6.88,8.00,2.51) will do the trick.
// Some city/years seem to stuck with a bad local optima with strange looking estimated distribution. Changing starting value to (0.95, -7.91, 11.8, -0.73) will do the trick. Among those include id 45, 2002, 2005, id 73, year 2001, 2006 and id 77, year 2002, 2006, id 25, year 2001.
// id 45 year 2006 and id 77 year 2001 didn't converge with multiple attempts with different starting values.

forv idval=1/111{
forv yearval=2001/2010{
use "./data/pm10data.dta", clear
qui keep if id==`idval' // city by city
qui keep if year==`yearval' // year by year
if _N>100{
scalar cff1=0.15
sum pmconc if pmconc<=float(cff1)
gen Fxnp=r(N)/_N  
gen cens=0

*****************
* CENSORED MLE *
*****************
scalar lftcff1=0.135   
scalar rtcff1=0.18
replace cens=1 if pmconc>float(lftcff1) & pmconc <=float(cff1)  
//gb2lfit_mod pmconc, cdf(pmparacdf2)  censvar(cens) 

scalar lna0=0.95
scalar lnb0=-7.91
scalar lnp0=11.8
scalar lnq0=-0.73
matrix intv = (lna0, lnb0, lnp0, lnq0)
matrix colnames intv = ln_a:_cons ln_b:_cons ln_p:_cons ln_q:_cons
gb2lfit_mod pmconc, cdf(pmparacdf2)  censvar(cens) from(intv) 

local ahat=e(ba)
local bhat=e(bb)
local phat=e(bp)
local qhat=e(bq)
local atx=cff1
gen Fxp=ibeta(`phat',`qhat', (`atx'/`bhat')^`ahat'/(1+(`atx'/`bhat')^`ahat'))
gen propFx=(Fxnp-Fxp)/Fxnp

local cityname=city[1] 
cumul pmconc, gen(pmecdf) 
gen cutoff=pmconc*0+cff1
line pmecdf pmconc if pmconc<=0.4, sort xlabel(0(0.1)0.4) lcolor(black) graphregion(color(white)) ///
|| line pmparacdf pmconc if pmconc<=0.4, lcolor(navy) sort ///
|| line pmecdf cutoff, lcolor(black) lpattern(dot) title("`cityname': `yearval'")  ///
name(`cityname'`yearval'_CMLEfit, replace) aspect(0.8) ///
legend(order(1 "Reported Data" 2 "Estimated Counterfactual Distribution") size(4) region(c(none)) bm(zero)) ytitle("") 
graph export "./outputgraphs/app/`cityname'`yearval'_CMLEfit.pdf", replace

keep in 1
if `idval'==1 & `yearval'==2003{
save "./data/PM10CMLEresults_Stata16.dta", replace
}
else {
append using "./data/PM10CMLEresults_Stata16.dta"
save "./data/PM10CMLEresults_Stata16.dta", replace
} 
}
}
}
